library(prism)
library(raster)
library(lubridate)
library(geosphere)
library(RColorBrewer)
library(png)
library(grid)
library(rasterVis)
library(rgdal)
library(scales)
library(viridis) 
library(ggthemes)
library(tidyverse)
library(lubridate)
library(DomoR)
library(choroplethr)

#get_prism_dailys(type="tmax", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="tmin", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="tmean", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="ppt", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)

retailCal <- read.csv("C:/Users/sewing/OneDrive - Do My Own/NOAA/4-4-5-calendar.csv",
                      stringsAsFactors = FALSE) %>% 
  mutate(Actual.Date = date(mdy(Actual.Date)),
         X2018.equiv = date(mdy(X2018.equiv))) %>%
  filter(Retail.Year > 2016,
         Retail.Year < 2019) %>%
  select(Actual.Date, Retail.Week, Retail.Year, Retail.DayofYear) %>%
  spread(Retail.Year, Actual.Date) %>%
  mutate(meanFile2017 = character(1),
         meanFile2018 = character(1),
         meanPath2017 = character(1),
         meanPath2018 = character(1)) %>%
  filter(!is.na(`2017`))

RetailWeeks2018 <- read.csv("C:/Users/sewing/OneDrive - Do My Own/NOAA/4-4-5-calendar.csv",
                      stringsAsFactors = FALSE) %>% 
  filter(Retail.Year == 2018) %>%
  select(Retail.Week, Actual.Date) %>%
  mutate(Actual.Date = mdy(Actual.Date)) %>%
  group_by(Retail.Week) %>%
  summarise(max_date = max(Actual.Date),
            min_date = min(Actual.Date)) %>%
  mutate(min_month = as.character(month(min_date, label = TRUE)),
         max_month = as.character(month(max_date, label = TRUE)),
         months = ifelse(min_month == max_month, min_month, paste0(min_month, "/", max_month))) %>%
  select(Retail.Week, months)

options(prism.path = "C:/Users/sewing/OneDrive - Do My Own/NOAA/PRISM_Files")
prismFileList <- ls_prism_data(absPath = TRUE)

stableMeans <- prismFileList[which(grepl("_tmean_stable_4kmD1_",
                                         prismFileList$files)),] %>%
  mutate(Date = ymd(str_sub(files, start = 26, end = -5)))

lastStableDate <- max(stableMeans$Date)

get_prism_dailys(type="tmean",
                 minDate = lastStableDate,
                 maxDate = Sys.Date(),
                 keepZip=F)

prismFileList <- ls_prism_data(absPath = TRUE)

stableMeans <- prismFileList[which(grepl("_tmean_stable_4kmD1_",
                                         prismFileList$files)),] %>%
  mutate(Date = ymd(str_sub(files, start = 26, end = -5)))

newLastStableDate <- max(stableMeans$Date)

provisionalMeans <- prismFileList[which(grepl("_tmean_provisional_4kmD1_",
                                              prismFileList$files)),] %>%
  mutate(Date = ymd(str_sub(files, start = -12, end = -5))) %>%
  filter(Date > ymd(newLastStableDate))

lastProvisionalDate <- max(provisionalMeans$Date)

earlyMeans <- prismFileList[which(grepl("_tmean_early_4kmD1_",
                                              prismFileList$files)),] %>%
  mutate(Date = ymd(str_sub(files, start = -12, end = -5))) %>%
  filter(Date > ymd(lastProvisionalDate))

lastEarlyDate <- max(earlyMeans$Date)

meansFiles <- stableMeans %>%
  rbind(provisionalMeans) %>%
  rbind(earlyMeans)

prismComps <- retailCal %>%
  mutate(meanFile2017 = meansFiles$files[which(meansFiles$Date == retailCal$`2017`)],
         meanFile2018 = ifelse(`2018` > lastEarlyDate,
                               "Future Date",
                               meansFiles$files[which(meansFiles$Date == retailCal$`2018`)]),
         meanPath2017 = meansFiles$abs_path[which(meansFiles$Date == retailCal$`2017`)],
         meanPath2018 = ifelse(`2018` > lastEarlyDate,
                               "Future Date",
                               meansFiles$abs_path[which(meansFiles$Date == retailCal$`2018`)]))
  
dailyRasters <- prismComps %>%
  mutate(`2017 Rasters` = list(length = 1),
         `2018 Rasters` = list(length = 1))

i <- 1

for(i in 1:dim(dailyRasters)[1]){
  dailyRasters$`2017 Rasters`[i] <- raster(dailyRasters$meanPath2017[i])
  i = i + 1
} 
 
i <- 1 

for(i in 1:table(dailyRasters$`2018`<lastEarlyDate)[2]){
dailyRasters$`2018 Rasters`[i] <- raster(dailyRasters$meanPath2018[i])
i = i + 1
}

anomCalc <- function(x, y) {
  return(x - y)
}

avgWeek <- function(a, b, c, d, e, f, g) {
  return((((a+b+c+d+e+f+g)/7)*(9/5))+32)
}
i <- 0
while(i < week(lastEarlyDate)-1){
  i = i+1
  weekOverlay17 <- overlay(dailyRasters$`2017 Rasters`[[(7*(i-1))+1]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+2]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+3]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+4]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+5]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+6]],
                          dailyRasters$`2017 Rasters`[[(7*(i-1))+7]],
                          fun = avgWeek)
  
  writeRaster(weekOverlay17,
              filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_2017_avg_temp"),
              overwrite=TRUE)
  
  weekOverlay18 <- overlay(dailyRasters$`2018 Rasters`[[(7*(i-1))+1]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+2]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+3]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+4]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+5]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+6]],
                          dailyRasters$`2018 Rasters`[[(7*(i-1))+7]],
                          fun = avgWeek)
  
  writeRaster(weekOverlay18,
              filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_2018_avg_temp"),
              overwrite=TRUE)
  
  weekYoYDiff <- overlay(weekOverlay18, weekOverlay17, fun = anomCalc)
  
  writeRaster(weekYoYDiff,
              filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_YoY_diff_avg_temp"),
              overwrite=TRUE)

}

YoYFiles <- dir(path = "C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/",
                pattern = "YoY_diff_avg_temp.grd",
                include.dirs = TRUE,
                full.names = TRUE)

us_map <- map_data("state")
i <- 1

while(i < week(Sys.Date())){
  r <- raster(paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_",
                        i,
                        "_YoY_diff_avg_temp.grd"))
  r_spdf <- as(r, "SpatialPixelsDataFrame")
  r_df <- as.data.frame(r_spdf)
  colnames(r_df) <- c("value", "x", "y")
  
  rp <- ggplot() +
    geom_path(
      data = us_map,
      aes(x = long,
          y = lat,
          group = group),
      color = "white",
      size = 0.1) +
    coord_equal() +
    theme_map() +
    theme(legend.position = "bottom") +
    theme(legend.key.width = unit(2, "cm")) +
    labs(
      x = NULL,
      y = NULL,
      title = paste0(RetailWeeks2018[i, 2], " - Week ", i, " YoY Diff in Temperature")
    ) +
    geom_raster(data = r_df, aes(x = x, y = y, fill = value)) + 
    scale_fill_distiller(palette = "RdYlBu", limits = c(-30, 30)) 
  
  fname <- paste0('C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Plots/Week_',
                  i,
                  '_2018_less_2017_Plot.png')
  ggsave(rp, file=fname)
  i = i+1
}